function [gbest,gbestval,fitcount] = MGWO_func(fhd,Dimension,Particle_Number,Max_Gen,lowerbound,upperbound,pos,bool_Plot,varargin)

rand('state',sum(100*clock));
me = Max_Gen;
ps = Particle_Number;
D = Dimension;

[pub, plb, vub, vlb] = SetBound(D,ps,lowerbound,upperbound);
% pos = VRmin+(VRmax-VRmin).*rand(ps,D);

vel = zeros(ps,D); % initialize the velocity of the particles
fitcount = ps;

fitness = feval(fhd,pos',varargin{:});
sort_fitness = [fitness' (1:ps)'];
sort_fitness = sortrows(sort_fitness, 1);
fitbest = sort_fitness(1,1);
fitworst = sort_fitness(ps,1);
gbestindex = sort_fitness(1,2);
gbest = pos(gbestindex,:);
gbestval = fitbest;

if(bool_Plot)
    figure(6)
    haxes = plot( 0 , 0 );
    XArray = [1];
    YArray = [gbestval];
    title('MGWO');
end

for iter = 2:me
    a=2-iter*((2)/me);
    Alpha_pos = pos(sort_fitness(1,2),:);
    Beta_pos = pos(sort_fitness(2,2),:);
    Delta_pos = pos(sort_fitness(3,2),:);
    
    for i = 1 : ps
        for j = 1 : D
            r1=rand(); % r1 is a random number in [0,1]
            r2=rand(); % r2 is a random number in [0,1]
            
            A1=2*a*r1-a; % Equation (3.3)
            C1=2*r2; % Equation (3.4)
            
            D_alpha=abs(C1*Alpha_pos(j)-pos(i,j)); % Equation (3.5)-part 1
            X1=Alpha_pos(j)-A1*D_alpha; % Equation (3.6)-part 1
                       
            r1=rand();
            r2=rand();
            
            A2=2*a*r1-a; % Equation (3.3)
            C2=2*r2; % Equation (3.4)
            
            D_beta=abs(C2*Beta_pos(j)-pos(i,j)); % Equation (3.5)-part 2
            X2=Beta_pos(j)-A2*D_beta; % Equation (3.6)-part 2       
            
            r1=rand();
            r2=rand(); 
            
            A3=2*a*r1-a; % Equation (3.3)
            C3=2*r2; % Equation (3.4)
            
            D_delta=abs(C3*Delta_pos(j)-pos(i,j)); % Equation (3.5)-part 3
            X3=Delta_pos(j)-A3*D_delta; % Equation (3.5)-part 3             
            
            vel(i,j)=(X1+X2+X3)/3;% Equation (3.7)
        end
    end
    pos = vel;
    pos = ((pos>=plb)&(pos<=pub)).*pos+(pos<plb).*plb+(pos>pub).*pub;
    
    fitness = feval(fhd,pos',varargin{:});
    sort_fitness = [fitness' (1:ps)'];
    sort_fitness = sortrows(sort_fitness, 1);
    fitbest = sort_fitness(1,1);
    fitworst = sort_fitness(ps,1);
    gbestindex = sort_fitness(1,2);
    gbest = pos(gbestindex,:);
    gbestval = fitbest;
    
    if(bool_Plot & mod(iter,10)== 0)
        XArray = [ XArray iter]; 
        YArray = [ YArray gbestval];
        set( haxes , 'XData' , XArray , 'YData' , YArray );
        drawnow
    end
end

end
